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Abstract. We present our latest results for simulation for merger of black hole (BH)- 
neutron star (NS) binaries in full general relativity which is performed preparing a 
quasicircular state as initial condition. The BH is modeled by a moving puncture with 
no spin and the NS by the T-law equation of state with T = 2 and corotating velocity 
field as a first step. The mass of the BH is chosen to be « 3.2M Q or 4.0M Q , and the 
rest-mass of the NS w 1.4M© with relatively large radius of the NS « 13-14 km. The 
NS is tidally disrupted near the innermost stable orbit but ~ 80-90% of the material 
is swallowed into the BH and resulting disk mass is not very large as ~ 0.3M Q even for 
small BH mass ~ 3.2M Q . The result indicates that the system of a BH and a massive 
disk of ~ Mq is not formed from nonspinning BH-NS binaries irrespective of BH mass, 
although a disk of mass ~ 0.1M© is a possible outcome for this relatively small BH 
mass range as ~ 3-4M©. Our results indicate that the merger of low-mass BH and NS 
may form a central engine of short-gamma-ray bursts. 



PACS numbers: 04.25.Dm, 04.30.-w, 04.40.Dg 
1. Introduction 

Merger of black hole (BH)-neutron star (NS) binaries is one of likely sources of kilo- meter 
size laserinterferometric gravitational wave detectors. Although such system has not 
been observed yet in contrast to NS-NS binaries, statistical studies based on the stellar 
evolution synthesis suggest that the merger will happen more than 10% as frequently as 
the merger of binary NSs [U l2j- Thus, the detection of such system will be achieved by 
laserinterferometers in near future. This motivates theoretical studies for the merger of 
BH-NS binaries. 

According to a study based on the tidal approximation (which is referred to as 
a study for configuration of a Newtonian star in circular orbits around a BH in its 
relativistic tidal field; e.g., 0, HI El El E] ) , the fate is classified into two cases, depending 
on the mass ratio q = M N g/M B H, where M B h and M N g denote the masses of BH and 
NS, respectively. For q <, q c , the NS of radius R will be swallowed into the BH horizon 
without tidal disruption before the orbit reaches the innermost stable circular orbit 
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(ISCO) [SI E]. On the other hand, for q > q c , NS may be tidally disrupted before 
plunging into BH. Here, the critical value of q c depends on the BH spin and equation of 
state (EOS) of NS. For the nonspinning case with stiff EOSs, q c « 0.3-0.35(i?/5M NS )~ 3/2 
and for the case that spin of a BH aligns with the orbital angular momentum, q c can be 
larger [6]. (Throughout this paper, we adopt the geometrical units c = G — 1.) 

The tidal disruption has been studied with great interest because of the following 
reasons, (i) Gravitational waves at tidal disruption will bring information about the 
NS radius since the tidal disruption limit depends sensitively on it |8J. The relation 
between the mass and the radius of NSs may be used for determining the EOS of high 
density matter [9]. (ii) Tidally disrupted NSs may form a massive disk of mass ~ 0.1- 
1M Q around BH if the tidal disruption occurs outside the ISCO. Systems consisting of 
a rotating BH surrounded by a massive, hot disk have been proposed as one of likely 
sources for the central engine of gamma-ray bursts (GRBs) with a short duration [10J, 
and hence, merger of low-mass BH and NS is a candidate. 

However, the scenario based on the tidal approximation studies may be incorrect 
since gravitational radiation reaction and self-gravitational effects of NS to the orbital 
motion are ignored. Radiation reaction shortens the time available for tidally disrupting 
NSs. The gravity of NS could increase the orbital radius of the ISCO and hence 
the critical value of the tidal disruption, q c , may be larger in reality. Miller [11] 
estimates these ignored effects and suggests that NSs of canonical mass and radius 
will be swallowed into BH without tidal disruption. Moreover, NSs are described by the 
Newtonian gravity in the tidal approximation. If we treat it in general relativity, the 
self-gravity is stronger and hence tidal disruption is less likely. 

Tidal disruption of NSs by a BH has been investigated in the Newtonian [12] and 
approximately general relativistic (GR) simulation [131 [7]. However, the criterion of the 
tidal disruption will depend on GR effects as mentioned above, and hence, a simulation 
in full general relativity is obviously required (see [2] for an effort). In [15], we present 
our first numerical results for fully GR simulation, performed by our new code which 
has been improved from previous one [El [17]; we handle an orbiting BH adopting the 
moving puncture method, which has been recently developed by two groups [18] (see 
also [19] for detailed calibration of this method). As the initial condition, we prepare 
a quasicircular state computed in a new formalism described in Sec. 2. Focusing 
particularly on whether NSs of realistic mass and radius is tidally disrupted to form 
a massive disk around nonspinning BHs, we illustrate that a disk with mass ~ M & is an 
unlikely outcome for plausible values of NS mass and radius and for BH mass greater 
than 3M & , although a disk of mass of O(0. 1M ) is possible. 

In this paper, we extend previous study; we perform a simulation for different BH 
mass from that in [15J to find the dependence of the disk mass on the BH mass. In 
addition, gravitational waves are derived. The paper is organized as follows. In Sec. 
2, we describe a formulation for computing quasicircular states in the moving-puncture 
framework. Section 3 presents some of numerical results for the quasicircular orbits. In 
Sec. 4, we report numerical results of simulations for the merger of BH-NS binaries. 
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Section 5 is devoted to a summary and discussion. 



2. Formalism for a quasicircular state 

Three groups have worked in computing quasicircular states of BH-NS binaries 
[201 I2T1 [22]. However, this field is still in an early stage in contrast to computation for 
NS-NS binaries (e.g., [23]). In [15] . we proposed a new method for computing accurate 
quasicircular states that can be used for numerical simulation in the moving puncture 
framework [2H[TH]. Here, we briefly describe the formulation. 

Even just before the merger, it is acceptable to assume that BH-NS binaries are 
in a quasicircular orbit since the time scale of gravitational radiation reaction is a few 
times longer than the orbital period. Thus, we assume the presence of a helical Killing 
vector around the mass center of the system, = (dt)^ + £1(8^)^, where the orbital 
angular velocity Q is constant. 

In the present work, we assume that the NS is corotating around the mass center 
of the system for simplicity. Irrotational velocity field is believed to be more realistic 
for BH-NS binaries [25]. The work for the irrotational case will be reported in the 
future paper [26]. The assumption of corotating velocity field in the helical symmetric 
spacetime yields the first integral of the Euler equation, h~ x u l = const, where h is 
specific enthalpy defined by 1 + e + P/p, and e, P, and p are specific internal energy, 
pressure, and rest-mass density, respectively. In the present work, we adopt the T-law 
EOS with r = 2; P = pe = Kp 2 with k an adiabatic constant, denotes the four 
velocity and u l its time component. Assumption of corotation implies = u 1 ^ and 
thus u r = u 6 = 0. 

For a solution of geometric variables of quasicircular orbits, we adopt the conformal 
flatness formalism for three-geometry. In this formalism, the solution is obtained by 
solving Hamiltonian and momentum constraint equations, and an equation for the time 
slicing condition which is derived from K k k = where is the extrinsic curvature 
and K k k its trace |27J. Using the conformal factor if), the rescaled tracefree extrinsic 
curvature = ip^K^ , and weighted lapse $ = aip where a is the lapse function, these 
equations are respectively written 

Aif) = -27rp H ^ 5 - \k^k;iT\ (1) 
A/, = 8vr J^ 6 , (2) 
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ij\p n + 2S) + —^A/A/ , (3) 



where A denotes the flat Laplacian, p H = ph^au 1 ) 2 — P, Ji = phui, and S = 
ph^au 1 ) 2 - 1] + 3P. 

We solve these equations in the framework of the puncture BH [2^| IT8| [28] . 
Assuming that the puncture is located at rp, we set if) and $ 

ijj = l + + and $ = 1 - — + 77, (4) 
2r B H r Bil 



Merger of black hole-neutron star binaries in full general relativity 4 

where M P and C are positive constants, and r B H = \ x bk\ ( x bh = xk ~ x p)- Then elliptic 
equations for functions <ft and 77 are derived. The constant Mp is arbitrarily given, while 
C is determined from the virial relation (e.g., [29|) 

/ d&dS* = -l d i ipdS i = 2vrM, (5) 

Jr— >oo Jr—*oo 

where M is the ADM mass. The mass center is determined from the condition that the 
dipole part of ip at spatial infinity is zero. In this method, the region with a < exists. 
However, this does not cause any pathology in the initial value problem. 
Equation (J2J) is rewritten setting 

4(= A^) = W id + W jti - 2 -5 l3 5 kl W Kl + Kf 3 , (6) 

where Kj- denotes the weighted extrinsic curvature associated with linear momentum 
of a puncture BH; 

K ij = ^~T~ ( n i P 3 + n J P i + ( n i n 3 ~ ^jO^fe) • ( 7 ) 
BH ^ ' 

Here, n k = n k = ^bh/^bh- P% denotes linear momentum of the BH, determined from 
the condition that the total linear momentum of system should be zero; 

Jii(j 6 d 3 x. (8) 

The RHS of Eq. (jSJ) denotes the total linear momentum of the companion NS. Then, 
the total angular momentum of the system is derived from 

J = J J^cPx + e^riS^Pt. (9) 

The elliptic equation for H / i (= W l ) is 

ISWi + l -d l d k W k = 8ttJ^ 6 . (10) 

Denoting W{ = 7B { — (x,i + B k ^x k ) where x an d B% are auxiliary functions [20], Eq. 
ffTUj) is decomposed as two linear elliptic equations 

ABi = irJii) 6 and Ax = -vrJixV 6 . (11) 

Computing BH-NS binaries in a quasicircular orbit requires to determine the shift 
vector even in the puncture framework. This is because has to be obtained [it is 
derived from u k = dki^ijj^iv 1 + j3 l ) where v % = ^l(f i \. The relation between f3 l and A\j is 
written 

3 W = j G 

Operating S^di, an elliptic equation is derived 

A/3* + \s tk d k d j f3 j = 2<9 i (aV>~ 6 )i* j + l&iraJjF*. (13) 

Here for A\j, we substitute the relation of Eq. (jBj) (not Eq. f[T2"j) ). As a result, no singular 
term appears in the RHS of Eq. ({TBI) and Eq. ({TBI) is solved in the same manner as that 
for W 4 . 



Sjkd^ + 5 lk d,(3 k - -6 {j d k p k = ^4, (12) 
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Table 1. Parameters of quasicircular states. Mass parameter of puncture, mass of BH, 
rest-mass of NS, mass, radius, and normalized mass of NS in isolation, total mass of 
the system, nondimensional angular momentum parameter, orbital period in units of 
M, and compactness of the system defined by C — (MQ) 2 / 3 . Mass of BH is computed 
from the area of the apparent horizon A as (A/I671-) 1 / 2 . Mass is shown in units of Mq. 





M P 


M BH 


M* 




R (km) 


M*/k 1 ' 2 


M 


J/M 2 


Po/M 


C 


A 


3.13 


3.21 


1.40 


1.30 


13.8 


0.147 


4.47 


0.729 


119 


0.141 


B 


3.13 


3.21 


1.40 


1.30 


13.8 


0.147 


4.47 


0.720 


110 


0.150 


C 


3.93 


4.01 


1.40 


1.30 


13.0 


0.151 


5.26 


0.645 


115 


0.144 



We have computed several models of quasicircular states and found that the relation 
between Q and J approximately agrees with the 3rd Post-Newtonian relation [21] . This 
makes us confirm that this approach is a fair way for preparing quasicircular states. 
We also found that in this method, the shift vector at r = rp automatically satisfies 
the condition (3^ = —Q within the error of a few %. This implies that the puncture is 
approximately guaranteed to be in a corotating orbit in the solution. 

3. Numerical results for quasicircular states 

BH-NS binaries in quasicircular orbits have been computed for a wide variety of models 
with q = M^/Mbh ~ 0.3-0.5 where M* denotes baryon rest-mass of the NS. In the 
present work, the compactness of spherical NSs with rest-mass M* is chosen to be 
~ 0.14-0.15. In the T-law EOSs, the mass and radius of NS are rescaled by changing 
the value of k: In the following we fix the unit by setting that M* = 1.4M . 

In computation, we focus only on the orbit of slightly outside of ISCO. In Table 
I, we show the quantities for selected quasicircular states with q « 0.4 and 0.33. For 
models A and B shown in Table I, the radius of NS in isolation is R ~ 13.8 km, the 
gravitational mass in isolation is M 0N s « 1.3OM , and MJk 1 ' 2 = 0.147. For model 
C, R « 13.0 km, M 0NS « 1.3M , and MJk 1 / 2 = 0.151. According to theories for NSs 
based on realistic nuclear EOSs [32], the radius of NS of Mns ~ 1-4M is 11-13 km. 
Thus the radius chosen here is slightly larger than that of realistic NSs and is more 
subject to tidal disruption. Model A is that used for the following numerical simulation 
and model B is very close to the tidal disruption limit of approximately the same mass 
as that of model A, showing that the model A has an orbit of slightly outside of the 
tidal disruption limit. This is also the case for model C. As we show in the next section, 
tidal disruption sets in after a small decrease of orbital separation for models A and C. 

The tidal approximation studies suggest that for q > q*, NSs could be tidally 
disrupted by a BH [6|. Here, in the tidal approximation, the critical value q* for T = 2 
and for nonspinning BHs is approximately given by 

q * = °H^- s ) — wr*—' (14) 
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and fl = is the angular velocity of the ISCO around nonspinning BHs. For 

models A and B, q* « 0.32, and hence, q > g*. According to the tidal approximation 
studies [HIE], such NS should be unstable against tidal disruption. Nevertheless, such 
equilibrium exists, proving that the tidal disruption limit in the framework of the tidal 
approximation does not give correct answer. Our studies indicate that the critical 
value q* is w 0.43( J R/5M NS )- 3 / 2 [(M BH fi)' 1 /6 3/2 ]; tidal disruption of NS is much less 
likely than in the prediction by the tidal approximation [5], [6]. For the typical NS 
of radius R ~ 5Mns and mass Mns ~ 1.4M , Mbh ^ 3.3M Q will be necessary for 
(Mbh^)~ 2 ^ 3 > 6; this implies that canonical NSs will not be tidally disrupted outside 
ISCO by most of nonspinning BHs of mass larger than ~ 3M . Tidal disruption occurs 
only for NSs of relatively large radius and only for orbits very close to ISCO. 

In this work, the criterion for the tidal disruption is investigated only for T = 2 
EOS and for compactness 0.14-0.15. The criterion is likely to depend on the stiffness of 
the EOS [6] since the structure of NS depends on it. The criterion should also depend 
strongly on the compactness of NS in general relativity which has a nonlinear nature. 
In the future paper, we plan to determine the criterion for a wide variety of EOSs and 
compactness of NS [26J. 



® = V + T-T^ ^ m , ■, 1 w » ( 15 ) 



4. Simulation for merger 

Even if tidal disruption of an NS occurs near ISCO, a massive disk may be formed around 
the companion BH. To investigate the outcome after merger and resulting gravitational 
waveforms, we perform numerical simulation adopting models A and C. 

For the simulation, we initially reset the lapse (i.e., $) since the relation a > 
should hold everywhere. In the present work, $ at t = is given by 

1 + 0.1X 4 
1 + Em=i X m + 1.1A 4 
where X = C/r^n- Then, a = only at puncture and otherwise a > 0. Furthermore, 
for r bh > C, the values of $ quickly approach to those of the quasicircular states. 

The numerical code for hydrodynamics is the same as that for performing merger 
of NS-NS binaries (a high-resolution central scheme) [331 [17] . On the other hand, we 
change equations for a, f3 l , and ip, and numerical scheme of handling the transport 
terms of evolution equations for geometries. For a and (3 % we solve 

(d t -[3%)\na = -2K k k , (16) 
d t (3 l = 0.75f j (F J +Atd t F J ), (17) 

where 7^- is the conformal three-metric and Fi = ^ 6^7^. At denotes the time step for 
the simulation and the second term in the RHS of Eq. (TTT|) is introduced for stabilization. 
The equation for the conformal factor is also changed to 

d t ^- G - d^- & (3 1 ) = (aK k k - 29#)^" 6 , (18) 

since ip diverges at the puncture [IB] . 
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In addition, we have improved numerical scheme for the transport term of geometric 
variables (d t — f3 l di)Q where Q is one of the geometric variables: First, we rewrite this 
term to d t Q — di{Q(3 l ) + Qdif3 l and then apply the same scheme as in computing the 
transport term of the hydrodynamic equations to the second term (3rd-order piece-wise 
parabolic interpolation scheme [IE]). We have found that for evolving BHs, such high- 
resolution scheme for the transport term in the geometric variables are crucial. This is 
probably because of the fact that near punctures, some of geometric variables steeply 
vary and so is the term (3 l diQ. For other terms in the Einstein's equation, we use the 
2nd-order finite differencing as in [T6J [17] . (Note that in the case of nonuniform grid, 
4-point finite differencing is adopted for Q,u since 3-point one is lst-order.) After we 
performed most of runs, we iterated some of computations with a 3rd-order scheme (5- 
point finite differencing for Q tii ) which is used in [HI [19]. We find that with such scheme, 
convergent results are obtained with a relatively large grid spacing. However, the results 
are quanlitatively unchanged and the extrapolated results (which are obtained in the 
limit of zero grid spacing; see below) are approximately identical. 

In the simulation, the cell-centered Cartesian, (x, y, z), grid is adopted to avoid the 
situation that the location of punctures (which always stay in the z = plane) coincides 
with the grid location. The equatorial plane symmetry is assumed and the grid size is 
(2N,2N, N) for x-y-z where N is a constant. Following [31], we adopt a nonuniform 
grid; in the present approach, a domain of (2N , 2N , N ) grid zone is covered with a 
uniform grid of the spacing Ax and outside the domain, the grid spacing is increased 
according to £ tanh[(z — iVo) / Az] Ax where % denotes the i-th grid point in each direction. 
iV , Az, and £ are constants. For model A, Ax/M P is chosen to be 1/8, 9/80, 1/10, 7/80, 
and 3/40, and for model C, it is 7/120, 8/120, 9/120, and 1/12. As shown in [18], such 
grid spacing can resolve moving punctures. For model A, (N, N , Ai, £, Ax/M p , L/X) 
is chosen to be (160,105,30,4.5,1/8,0.46), (200,105,30,4.5,1/8,0.78), (200,105,30,5,1/8, 
0.83), (220,125,30,5,9/80, 0.78), (220,125,30,6,1/10, 0.78), (220,140,30,7,7/80, 0.65), and 
(220,150,9,3/40, 0.59). Here, L and A denote the location of the outer boundaries along 
each axis and the wavelength of gravitational waves at t — 0. For model A with 
Ax = Mp/8 and N = 160, we chose other values of N and Ai, and found that results 
depend weakly on them as well as on L as far as L ^ A/2. For model C, the chosen 
parameters are (200, 120, 30, 6, 1/12, 0.53), (220, 125, 30, 6, 9/120, 0.58), (220, 125, 
30, 7, 8/120, 0.58), and (220, 140, 30, 8, 7/120, 0.48). 

For a test, we performed simulations for merger of two nonsp inning BHs adopting 
the same initial condition as used in [18] . We focused particularly on the merger time, 
which is referred to as the time at which a common apparent horizon is first formed 
and found that it varies with improving grid resolution approximately at 1st order. 
The likely reason is that geometric variables vary steeply around the BH where they 
are evolved with lst-order accuracy in our scheme, although other regions are resolved 
with 2nd-order accuracy. By the extrapolation, an exact merger time is estimated to 
be ~ 19M. This result agrees approximately with those of [H]. This indicates that 
our code can follow moving punctures (see the left panel of Fig. [1] for evolution of the 
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Figure 1. Numerical results for merger of binary BHs performed with the initial 
data of [18]. Left: The open thick circles denote evolution of the location of the 
apparent horizon of two BHs for t/M = 0, 5, 10, 15, and 20 where M is the total ADM 
mass at t = 0. The open wide circle at t = 20M is the common apparent horizon. 
The small solid circles denote approximate location of the maximum of the conformal 
factor. The grid setting is (N, No, Ai, £, Ax/M p , L/X) = (200,105,30,7,1/16,1.20) 
for this result (note that A 37M). Right: + and x modes of gravitational 
waveforms extracted r « A for (N, No, Ai, £, Ax/M p , L/X) = (200,105,30,9,0.05,1.20) 
(solid curves), (200,105,30,7,0.064,1.23) (dashed curves), and (200,105,30,5,0.08,1.23) 
(dotted curves). 



location of apparent horizons for the initial condition of [18J). 

Gravitational waves are also computed. In the right panel of Fig. [T] we display 
the results for three different grid resolutions. The figure shows that after merger, 
gravitational waves are determined by the quasinormal mode ringing and that the 
waveforms depend very weakly for chosen grid resolutions. The wavelength of the 
quasinormal mode is ~ 11-12M, which agrees with the result of |18j . 

Next, we present the results for merger of BH and NS for model A. Figure [2] 
shows evolution of contour curves for p and velocity vectors for v % in the equatorial 
plane together with the location of apparent horizons at selected time slices for 
(N, N , Ai, £, Ax/M p , L/X) =(220,150,9,3/40, 0.59). Due to gravitational radiation 
reaction, the orbital radius decreases and then the NS is elongated (2nd panel). Because 
of the elongation, the quadrupole moment of the NS is amplified and the attractive force 
between two objects is strengthen [35]. This effect accelerates an inward motion and, 
consequently, the NS starts plunging to the BH at t ~ 90M. Soon after this time, the 
NS is tidally disrupted; but the tidal disruption occurs near the ISCO and hence the 
material in the inner part is quickly swallowed into the BH (3rd and 4th panels). On 
the other hand, because of the outward angular momentum transfer, the material in the 
outer part of the NS forms a disk with the maximum density ~ 10 12 g/cm 3 (5th and 
6th panels). 




Figure 2. Snapshots of the density contour curves for p in the equatorial plane 
for model A with (N, N Q ,Ai, £, Ax/M p , L/X) =(220,150,9,3/40, 0.59). The solid 
contour curves are drawn for p = (1 + 2i) x 10 14 g/cm 3 (i = 1,2,3) and for 
10 u-o.5i g / cm 3 (j = i ^ 8). The maximum density at t = is w 7.2 x 10 14 g/cm 3 . 
The blue, cyan, magenta, and green curves denote 10 14 , 10 12 , 10 11 , and 10 10 g/cm 3 , 
respectively. Vectors indicate the velocity field (v x ,v v ), and the scale is shown in the 
upper right-hand corner. The thick (red) circles are apparent horizons. Time is shown 
in units of total mass of the system M. 



Mass of the disk is not large. Figure G](a) shows evolution of baryon rest-mass 
located outside apparent horizons M r>rAu . We find that ~ 80% and ~ 90% of the mass 
is swallowed into the BH in t ~ 130M ~ 2 ms for models A and C, respectively, for the 
best grid-resolution simulations. The swallowing continues after this time. Note that 
since M r>rAn is rest-mass outside apparent horizons, the disk mass which should be 
defined for the mass located outside an ISCO around the formed BH is slightly smaller. 
We follow the evolution of the rest-mass of material located for r > 3M and r > 4.5M, 
where r = 3M-4.5M are approximate locations of the ISCO around the formed BH, 
and find that their values are smaller than M r>rAH by ~ 10% and 20%, respectively. 
Thus, the disk mass would be 0.8-0.9M r>rAH in reality. 

The value of M r>rAU depends systematically on Ax; we find that the results for 
good resolutions {Ax/M p = 3/40, 7/80, and 1/10 for model A and 7/120, 8/120, and 
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Figure 3. (a) Results for model A for various grid resolutions. Upper panel shows 
evolution of baryon rest-mass located outside the apparent horizon. The plots for 
Aa; = Mp / 8 almost coincide and show that the results depend weakly on the values 
of L and £. On the other hand, the results depend systematically on Aa: (see text). 
Lower panel shows evolution of area of apparent horizon in units of 16irM 2 . (b) The 
same as (a) but for model C. 



9/120 for model C) at late times approximately obey a relation of convergence, i.e., 
M r>rAH (t) = a(t) + b(t)Ax n where a(t) and b(t) are functions of time. The order of 
convergence, denoted by n, is between 1st and 2nd order (i.e., 1 < n < 2). For model A, 
least-square fitting gives a{t) at t = 180M as ~ 0.19M* if we set n = 2 and as 0.28M* 
for n = 1 (see the solid circle in Fig. [3]^a)). For model C, they are 0.11M* and 0.15M* 
(the solid circle in Fig. G^b)), respectively. Thus, the true result should be between 
0.19M* and 0.28M* for model A and between 0.11M* and 0.15M* for model C. 

It is reasonable that the disk mass for model C is smaller than for model A since 
the mass ratio of NS to BH is smaller for model C while the compactness of NS is nearly 
identical. A remarkable point is that with a small increase of the BH mass from 3.2M 
to 4.0M Q , the disk mass decreases by a factor of 2. This suggests that the disk mass 
depends sensitively on the BH mass. 

Note that the adopted NS has corotating velocity field initially and furthermore 
its radius is larger than that of canonical NSs. For irrotational velocity field with a 
realistic value of radius, the disk mass would be smaller than this value. It is reasonable 
to consider that our results for the disk mass provide an upper limit of the disk mass 
for given mass and spin of BH. Hence it is unlikely that a massive disk with ~ M & is 
formed after merger of nonspinning BH of mass M > 3M Q and canonical NS of mass 
m 1AM & and radius ~ 11-13 km, although a disk of mass of ~ 0.2-0.3M Q may be 
formed for small BH mass. This value of the disk mass is large enough to explain short 
GRBs of relatively low total energy ~ 10 49 ergs (e.g., [36J for an estimate). 

Lower panels of Fig. [3] show the evolution of the area of apparent horizons in units 
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of l&irM 2 . These illustrate that the area of BHs quickly increases swallowing material, 
and then, the area settles down to an approximate constant. We, again, evaluate the 
true final area by extrapolation for results with different grid resolutions. For model 
A, the area in units of I&ttM 2 at t — 180M is 0.74 and 0.80 in the assumption of the 
1st (n = 1) and 2nd-order (n = 2) convergences, respectively (see the solid circles in 
Fig. E(a)). For model C, A/lQirM 2 at t = 200M is 0.84 and 0.88 for n = 1 and n = 2, 
respectively (see the solid circles in Fig. [3(b)). From these values the spin parameter of 
formed BHs, a, is approximately derived from 



A 1 + VT- a 



167iM 2 m ~ 2 



(19) 



where Mem denotes the mass of the formed BH. To approximately estimate it, we 
simply use 

M BHf = M — M r>rAH - E GW} (20) 

where £"gw is radiated energy by gravitational waves. We find that -Egw is about 1% 
of M and simply use relation £"gw — 0.01M. Then we obtain a = 0.57 and 0.52 for 
n = 1 and 2 for model A and a = 0.52 and 0.42 for n = 1 and n = 2 for model C. Thus, 
spinning BHs of moderate rotation are outcomes. 

The spin parameter of the formed BHs is much smaller than the initial value 
of J/M 2 of the system. One of the reasons is that the disk has large angular 
momentum approximately written as 3MBHf^disk where Mdi S k denotes the disk mass 
~ 0.8-0.9M r>rAH Here, the factor 3M B Hf denotes a value of typical specific angular 
momentum around the formed BH. Denoting the initial angular momentum by a^M 2 
where a 0.73 and 0.65 for models A and C (see Table I), the fraction of the angular 
momentum that the disk has is ~ 3oq 1 M^s^M-Qm / M 2 . Thus, for model A, the fraction 
is ~ 20-30% and for model C, it is 10-15%. In addition, gravitational waves carry away 
the angular momentum by ~ 10% of a M 2 . Thus, the angular momentum of the formed 
BH should be smaller than the initial total angular momentum of the system by 30-40% 
for model A and by 20-25% for model C. Therefore, the values for a derived above are 
reasonable magnitudes. 

In Fig. HI gravitational waveforms for model C is shown. Gravitational waves 
are extracted from the metric near the outer boundaries using a gauge-invariant wave 
extraction method (see [TOJ [33] for details). From the values of R + and R x , the maximum 
amplitude of gravitational waves at a distance D is evaluated 



Here, the maximum amplitude can be observed if the observer is located along the z- 
axis. Figure H] implies that the maximum amplitude at a distance of D = 100 Mpc is 
w 5 x 10~ 22 since M = 5.26M . 

For t TCt <, 120M, inspiral waveforms are seen: Amplitude increases and 
characteristic wavelength decreases with time. The wavelength at the final phase of 
the inspiral is ~ 25M indicating that the orbital period of the ISCO is ~ 50M (i.e., 
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the angular velocity is ~ 0.12M -1 ). This is in approximate agreement with the 3rd 
post-Newtonian results |31j. 

For 120M <, tret ^ 150M, ring-down waveforms are seen. The characteristic 
wavelength is ~ 15M^nh which is in approximate agreement with the wavelength of 
the quasinormal mode. 

The characteristic feature of gravitational waves after the merger sets in is that 
the amplitude damps quickly even in the formation of massive disk. This is due to the 
fact that the degree of nonaxisymmetry of the disk decreases in a very short time scale 
(~ 20-30M). This indicates that in the frequency domain, the amplitude of Fourier 
power spectrum steeply decreases in the high-frequency region. 

Figure [5] shows the evolution of averaged violation of the Hamiltonian constraint. 
For the average, rest-mass density is used as a weight (see [16] for definition) and the 
integral is performed for the region outside apparent horizons. The figure shows that the 
Hamiltonian constraint converges approximately at 2nd order. This result is consistent 
with the fact that the region except for the vicinity of BH is followed with 2nd-order 
accuracy. 

5. Discussion 

In this paper, we have presented our latest numerical results of fully GR simulation for 
merger of BH-NS binary, focusing on the case that the BH is not spinning initially and 
the mass ratio q is fairly large as 0.3-0.4. It is found that even with such high values 
of q, the NS is tidally disrupted only for the orbit very close to ISCO and 80-90% of 
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Figure 5. Evolution of averaged violation of the Hamiltonian constraint (a) for model 
A and (b) for C. 

the mass element is quickly swallowed into the BH without forming massive disks. The 
results do not agree quantitatively with the prediction by the tidal approximation study 
The reasons are: (1) In the tidal approximation, one describes NSs by the Newtonian 
gravity. In general relativity, gravity is stronger and the tidal disruption is less likely. (2) 
The time scale for angular momentum transfer during tidal disruption near the ISCO is 
nearly as long as the plunging time scale determined by gravitational radiation reaction 
and attractive force between two objects. Hence before the tidal disruption completes, 
most of the material is swallowed. 

If the BH has a large spin, the final fate may be largely changed because of the 
presence of spin-orbit repulsive force. This force can weaken the attractive force between 
BH and NS and slow down the orbital velocity, resulting in smaller gravitational wave 
luminosity and longer radiation reaction time scale [3?]. This effect may help massive 
disk formation. The study of spinning BH binaries is one of the next issues. The fate 
will also depend on EOS of NS [6] and mass of BH as well as on the velocity field of 
NS. Simulations with various EOSs and BH mass and with irrotational velocity field are 
also the next issues. 

In this paper, a small number of the results for quasicircular states of BH-NS 
binaries are presented. Currently we are working in the computation of quasicircular 
states for a wide variety of masses of BH and NS in the framework described in Sec. 2. 
The numerical results will be reported in the near future [26] . 
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